function [a_q_solve] = a_q_calib_fn(inputs_struc,fns_struc)

% This function calibrates the a_q's to match the targeted share of
% employment in different quality levels. Given the wage primium that will
% be calibarted and the targeted M_q (comes from target for average size),
% we can figure out how much will have to be demanded by the consumer in
% order to generate the required employment in each quality level (i.e. to match the supply implied by employment targets). 
% A good guess for the a_q's can be contructed. This is then used as a
% strating point to solve for the exact a_q's.
% One of the a_q's can be normalized to 0.
% The other are calibrated to match the ratio of supplies/demands of the
% other q's relative to the worst quality q.

w_U = inputs_struc.w_U; w_S_target = inputs_struc.w_S_target; theta_q_i = inputs_struc.theta_q_i;
A_q_i = inputs_struc.A_q_i; sigma_us = inputs_struc.sigma_us; p_z = inputs_struc.p_z;
z_input_q_i  = inputs_struc.z_input_q_i; sigma = inputs_struc.sigma;
sh_emp_q_target = inputs_struc.sh_emp_q_target; entry_sh_L = inputs_struc.entry_sh_L;
N = inputs_struc.N; size_q_target = inputs_struc.size_q_target;
pdf_A_q_i  = inputs_struc.pdf_A_q_i; eta = inputs_struc.eta;
prod_grid_size=inputs_struc.prod_grid_size; theta_q = inputs_struc.theta_q;
q = inputs_struc.q; h = inputs_struc.h; pi_target = inputs_struc.pi_target;
A_S = inputs_struc.A_S;

MC_L_fn = fns_struc.MC_L_fn; L_S_fn = fns_struc.L_S_fn;

% Computing all prices etc. Note that the A_q_i is scaled so at to
% attain the desired price ratios across qualities. This is assuming we hit
% the wage and other targets in equi based on calibration of other
% parameters.
MC_L_q_i = MC_L_fn(w_U,w_S_target,theta_q_i,A_q_i,sigma_us,A_S);
p_q_i = (MC_L_q_i + p_z*z_input_q_i)*sigma/(sigma-1);
P_q_norm = (sum((p_q_i .^(1-sigma)).*pdf_A_q_i ,2)).^(1/(1-sigma));
M_q_target = sh_emp_q_target.*(1-entry_sh_L).*N./size_q_target;    % The mass of firms we need in equilibrium to match the targeted size of each q
P_q = (M_q_target.^(eta - 1/(sigma-1))).*(sum((p_q_i .^(1-sigma)).*pdf_A_q_i ,2)).^(1/(1-sigma));

% Some useful constants
L_q_constant = ((w_U/w_S_target*(1-theta_q)./theta_q).^sigma_us)*(A_S^(sigma_us-1));

% Figuring out what the supply of each quality will be given the entry
% and the share of employment in each q. See the file
% "FGH_CES_lnC_HetrProd_LoVParameter_calib.lyx" for more details
% (section titled "Utility parameters a_q's given q's")
L_U_q = sh_emp_q_target.*(1-entry_sh_L).*N./(1 + L_q_constant);     %How much unskilled labor will be employed by each quality
l_j_to_i_ratio = (A_q_i./repmat(A_q_i(:,1),1,prod_grid_size)).^(sigma-1);       % Ratio of unskilled labor for two different productivities within a q
l_j_q_U = L_U_q./M_q_target./(sum(l_j_to_i_ratio.*pdf_A_q_i,2));   % amount of unskilled labor employed by the lowest productivity of q, given total labor employed in that q.
l_U_q_i = repmat(l_j_q_U,1,prod_grid_size).*l_j_to_i_ratio; % Amount of unskilled labor employed by all productivities of all q's
l_S_q_i = L_S_fn(l_U_q_i,theta_q_i,w_U,w_S_target,sigma_us,A_S); % Amount of skilled labor employed by all productivities of all q's
x_q_i = A_q_i.*l_U_q_i.*(theta_q_i + (1-theta_q_i).*(((1-theta_q_i)./theta_q_i*w_U/w_S_target*A_S).^(sigma_us-1))).^(sigma_us/(sigma_us-1));    % Amount produced bby each q,A
Y_q_supply = 1./(M_q_target.^eta).*(M_q_target.*sum((x_q_i.^((sigma-1)/sigma)).*pdf_A_q_i,2)).^(sigma/(sigma-1)); % Amount produced by each final quality producer

% Computing a good guess for a_q's. This is based on the equation labeled
% 'a_q guess' in the file "FGH_CES_lnC_HetrProd_LoVParameter_calib.lyx".
mean_q = mean(q);
w_U = w_U + pi_target;
w_S_target = w_S_target + pi_target;
Y_ratio_term = (1-h)*w_U./P_q + h*((w_U/w_S_target)^mean_q)*((w_S_target/w_U).^q).*(w_S_target./P_q);
exp_a_q = Y_q_supply(1,1)./Y_q_supply*(P_q(1,1)^q(1,1)).*(P_q.^(-q)).*Y_ratio_term./(Y_ratio_term(1,1));
a_q_guess = -[log(exp_a_q)];

% Structure used for inputs into the fsolve for a_q's
inputs_aQ_struc.Y_q_supply = Y_q_supply;
inputs_aQ_struc.q =  q;
inputs_aQ_struc.w_U = w_U;
inputs_aQ_struc.w_S_target = w_S_target;
inputs_aQ_struc.h = h;
inputs_aQ_struc.P_q = P_q;

% Function handle which gets inputed into the fsolve.
a_q_zero_handle = @(a_q_norm)a_q_zero(a_q_norm,inputs_aQ_struc);

% Solves for the a_q's which match the share of employment in each q.
% options_a = optimset('Display','final','Jacobian','on','DerivativeCheck','on');
% options_a = optimset('Display','final','Jacobian','on');
% options_a = optimset('Display','final');
options_a = optimset('Display','off');
a_q_solve = fsolve(a_q_zero_handle,a_q_guess(2:end,1),options_a);



end


function [y,J] = a_q_zero(a_q_norm,inputs_aQ_struc)

% Takes a set of a_q's as given and returns 'y' which needs to be set to 0.
% For a given set of a_q's and q's, it finds the demand for each quality,
% and checkes to see how much is the divergence from the supply needed to match employment targets.
% It is used to calibrate the a_q's given q as the a_q's are set so that
% the demand equals the targeted supply.
% The file     % "FGH_CES_lnC_HetrProd_LoVParameter_calib.lyx" for more details
% (section titled "Utility parameters a_q's given q's")

Y_q_supply 	= inputs_aQ_struc.Y_q_supply;
q 			= inputs_aQ_struc. q;
w_U 		= inputs_aQ_struc.w_U;
w_S_target 		= inputs_aQ_struc.w_S_target;
h 			= inputs_aQ_struc.h;
P_q 		= inputs_aQ_struc.P_q;

a_q = [0;a_q_norm];
Y_q_supply_norm = Y_q_supply(2:end,1);  % Amount produced 
q_norm = q(2:end,1);
P_q_norm = P_q(2:end,1);

g = sum(exp(a_q).*((w_U./P_q).^q))/sum(exp(a_q).*((w_S_target./P_q).^q));

y_1 = log(Y_q_supply(1,1)./Y_q_supply_norm.*(w_U.^(q_norm - q(1,1))).*(P_q_norm.^(-q_norm))*(P_q(1,1).^q(1,1)));
y_2 = a_q_norm;
y_3 = log((1-h)*w_U./P_q_norm + h*((w_S_target/w_U).^q_norm).*g.*(w_S_target./P_q_norm));
y_4 = log((1-h)*w_U./P_q(1,1) + h*((w_S_target/w_U).^q(1,1)).*g.*(w_S_target./P_q(1,1)));

% This is the log of the demand for a given quality relative to base
% quality. This needs to be set to 0 basically. It is what is labelled 'a_q
% zero eq' in the file "FGH_CES_lnC_HetrProd_LoVParameter_calib.lyx". These
% are the f() equation which need to be set to 0 basically.
y = y_1 + y_2 + y_3 - y_4;

% This the analytical Jacobian to the f() equations. See
% "FGH_CES_lnC_HetrProd_LoVParameter_calib.lyx"
d_g_by_d_a = exp(a_q_norm).*((1./P_q_norm).^q_norm).*(w_U.^q_norm - g.*(w_S_target.^q_norm))./(sum(exp(a_q).*((w_S_target./P_q).^q)));
d_f_by_d_i_pt1      = h*((w_S_target/w_U).^(q_norm)).*(w_S_target./P_q_norm)./((1-h)*w_U./P_q_norm + h*((w_S_target/w_U).^(q_norm)).*(w_S_target./P_q_norm).*g);
d_f_by_d_i_pt1_norm = h*((w_S_target/w_U).^(q(1,1))).*(w_S_target./P_q(1,1))./((1-h)*w_U./P_q(1,1) + h*((w_S_target/w_U).^(q(1,1))).*(w_S_target./P_q(1,1)).*g);

J = eye(size(q,1)-1) + d_f_by_d_i_pt1*d_g_by_d_a' - d_f_by_d_i_pt1_norm*repmat(d_g_by_d_a',size(q,1)-1,1);


end



